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Abstract 



Let n be a positive integer and to be a positive even integer. Let 
A be an ttt,*'' order n-dimensional real weakly symmetric tensor and B 
be a real weakly symmetric positive definite tensor of the same size. 
r^ I A e IR is called a S^- eigenvalue of A if Ax"'^~^ — ASa;™^^ for some 

"jrt . X € ]R"\{0}. In this paper, we introduce two unconstrained optimiza- 

tion problems and obtain some variational characterizations for the 
minimum and maximum K^ -eigenvalues of A. These unconstrained 
optimization problems can be solved using some powerful optimiza- 
'^ . tion algorithms, such as the BFGS method. We provide some numer- 

ical results to illustrate the effectiveness of this approach for finding 
a Z-eigenvalue and for determining the positive semidefiniteness of an 
even order symmetric tensor. 
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1 Introduction 

Since the pioneering works of Qi [17J and Lim [T3], the tensor eigenprob- 
lem has become an important part of numerical multilinear algebra. In this 
paper, we consider the real eigenvalue problems for even order real sym- 
metric tensors. Eigenvalues of symmetric tensors have found applications in 
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several areas, including automatic control, statistical data analysis, higher 
order diffusion tensor imaging, and image authenticity verification, etc., see 
for example, HH HSl |20l [HI [22] . 

Throughout this paper, we assume that IR is the real field, n is a positive 
integer, and m is a positive even integer. An m^'^-order n-dimensional real 
tensor 



)nxnx-xn 



is called symmetric if its entries are invariant under any permutations of their 
indices [I0l[l7]. A tensor A is called positive definite (positive semidefinite) 
if the multilinear form 

n 
■A.X / ^ ■^iii2---im'^il'^i2' ' ' -^im 

is positive (nonnegative) for all x G ]R"\{0}. The notation Ax^~^ denotes 
the vector in ]R° whose i entry is 

n 
[A.X ji = / , -^ii2---im-^i2 ' ' ' ■'^im- 

i2,---,im = l 

Following [6] , ^ is called weakly symmetric if the gradient of Ax"^ 

V (^x™) = mAx'^-^ 

for all X G H". If A is symmetric, then it is weakly symmetric [6]. 

Various definitions of real eigenpairs for tensors have been introduced 
in the literature, including H-eigenvalues [TF], Z-eigenvalues [IT], and D- 
eigenvalues [20]. We use the following generalized eigenvalue definition, 
which includes the H-, Z-, and D-eigenvalues as special cases. 

DEFINITION 1 Let A and B he m -order n-dimensional real weakly 
symmetric tensors. Assume further that B is positive definite. If there exist 
a scalar A G H and a nonzero vector x G JR^ such that 

Ax"^-^ = XBx'^-\ (1.1) 

then A is called a Br ^eigenvalue of A and x a Br~eigenvector with respect 
to A. We denote the BrSpectrum of A by 

f'"Br(-^) = {A : A is a ^Br— eigenvalue of A}. 



REMARK 1 This definition was first introduced by Chang, Pearson, and 
Zhang [6] in a somewhat more general setting. 

• li B = I = {Sij^i2---im)^ the unit tensor, then B is weakly symmetric 
positive definite. Moreover, Bx"^ = \\x\\^ = x^ + X2^ + • • • + x^, 
Bx'^~^ = m[3;^"\x™~\--- ,x^"^]'^, and the ^^-eigenvalues are H- 
eigenvalues. 

• li B = I™ , the tensor product of m/2 copies of the unit matrix 
In £ K"^'^, then B is weakly symmetric positive definite. Moreover, 
Bx"^ = (x^x)™/2 = {xl + xl + --- + x2)"^/2, Bx-^-^ = m{x'^x)^x, 
and the ^S^-eigenvalues are Z-eigenvalues. 

• li B = D^'"^ , the tensor product of m/2 copies of the symmetric pos- 
itive definite matrix D G H'^^'^, then B is weakly symmetric positive 
definite. Moreover, Bx"^ = {x^Dx)"^^"^, Bx"^~^ = m{x^ Dx)^2- Dx, 
and the i3r.-eigenvalues are D-eigenvalues. 

Calculation of all eigenvalues of a high order (m > 2) tensor is difficult, 
unless m and n are small [T7]. In certain circumstances, however, one only 
needs to compute the largest or smallest eigenvalue of a tensor. For instance, 
the smallest H-eigenvalue or Z-eigenvalue of an even order symmetric tensor 
A can be used to determine the positive definiteness/semidefiniteness of A 
|17j . For a nonnegative tensor, the Perron- Frobenius theory asserts that its 
largest H-eigenvalue is its spectral radius [H |8] . 

Recently, Kolda and Mayo [11] have extended the high order power 
method for symmetric tensor eigenproblems of Kofidis and Regalia [lOj by 
introducing a shift parameter a to compute Z-eigenvalues of symmetric ten- 
sors. With a suitable choice of a, the resulting method, SSHOPM, converges 
to a Z-eigenvalue of the tensor when applied to a symmetric tensor. The 
found Z-eigenvalue is not necessarily the largest or smallest Z-eigenvalue. 
The rate of convergence of the SSHOPM method is linear [11]. 

An alternative approach for computing the eigenvalues of a symmetric 
tensor is to solve the constrained optimization problem J13j 



or 



mmAx"" s.t. Bx"" = 1, (1.2) 



max Ax"" s.t. Bx"" = 1. (1.3) 



The Karush-Kuhn- Tucker points of Problem (jl.2p or (jl.Sp give iSf-eigenvalues 
and i3j.-eigenvectors of A. If we are interested in obtaining one eigenvalue. 



then these problems can be solved using a local constrained optimization 
solver [15] . Note that in each problem, the objective function and the con- 
straint function are both polynomials. Therefore, a global polynomial opti- 
mization method can be used, if we are interested in finding the largest or 
smallest iS^-eigenvalue. 

A more attractive approach for computing eigenvalues of even order sym- 
metric tensors, however, is to use unconstrained optimization. This is mo- 
tivated by the works of Auchmuty [lil2j, in which he proposed some uncon- 
strained variational principles for generalized symmetric matrix eigenvalue 
problems. In particular, he [2] considered the unconstrained optimization 
problems 

min gi{x) = -{x^ Bxf + -x^ Ax, (1.4) 

xgR" 4 2 

and 

min qoix) = —{x Bx) Ax x, (1-5) 

where A G IR'^^'^ is a symmetric matrix and B G IR'^^'^ is a symmetric pos- 
itive definite matrix. He proved that Problem (jl.4p can be used to find the 
smallest generalized i?-eigenvalue of A and Problem (jl.Sp can be used to 
find the largest generalized i?-eigenvalue of A. In this paper, we will ex- 
tend Auchmuty's unconstrained variational principles for symmetric matrix 
eigenproblems [2] to even order weakly symmetric tensors. 

The rest of this paper is organized as follows. In Section 2, we introduce 
some preliminary results that will be used to establish the main results 
in Section 3. In Section 3, we introduce two unconstrained optimization 
problems and obtain some variational characterizations for the minimum 
and maximum iS^-eigenvalues of A. In Section 4, we give some numerical 
results. Some final remarks are given in Section 5. 

2 Preliminaries 

We start with the existence of i3r--eigenvalues of A. The existence of H- 
eigenvalues and Z-eigenvalues of an even order symmetric tensor was first 
studied by Qi [ITj. In [6j, Chang, Pearson, and Zhang proved the existence 
of at least n ;Sr.-eigenvalues when A is weakly symmetric and B is weakly 
symmetric positive definite, which is summarized in the following 

THEOREM 1 (16]) Assume that A and B are m^^-order n-dimensional 
real weakly symmetric tensors and B is positive definite. Then A has at 
least n Br -eigenvalues, with n distinct pairs of Br~eigenvectors. 



In |17j . Qi proved the existence of the maximum and minimum H- 
eigenvalues and Z-eigenvalues. Using a similar argument, we can prove 

THEOREM 2 Assume that A and B are m -order n-dimensional real 
weakly symmetric tensors and B is positive definite. Then cs^^A) is not 
empty. Furthermore, there exist Amin G crB^{A) and Amax S o-is^{A) .such 
that 

-OO < Amin < A < Amax < OO, V A G (Tl3r{A). 

Proof : Since B is positive definite, the set {x £ ]R° : fix™ = 1} is compact 
[6]. We also notice that function Ax"^ is continuous. Thus, the constrained 
optimization problem ()1.2p has a global minimizer x and the constrained 
optimization problem (jl.3p has a global maximizer x. 

At the global minimizer x of problem (ll.2p , there is a scalar A such that 
the KKT conditions 

mAx""-^ = XmBxT-^ (2.1) 

hold. Clearly, A G ctq^^A). The inner product of ()2.ip with x gives 

Ax^ = XBx"" = A. 

Since x is a global minimizer of problem (jl.2p , 

A<A, \/XeaBr{A). 

Therefore, we can set Amin = A- Similarly, we can establish the existence of 
Amax by using the global maximizer x. □ 

We next consider a property of weakly symmetric positive definite ten- 
sors, which generalizes a similar property for symmetric positive definite 
matrices. 

THEOREM 3 Assume that B is an m^^-order n-dimensional weakly sym- 
metric positive definite tensor. Let fj, > be the smallest H-eigenvalue of B. 
Then 

Bx"" > ;u||xC \/x e ]R°, (2.2) 

where \\x\\m is the m-norm of x. 

Proof : When a; = 0, (j2.2p obviously holds. According to Theorem [21 /x is 
the global minimum value of 

minBx'", s.t. ||x||;^ = 1. 



For any x £ ]R'^\{0}, we have 



^(^) >M. 



This imphes (j2:2)) . □ 

Finahy we recah that a continuous function / : IR"^ — )• ]R is coercive if 

hm f{x) = +00. 

|x||— >oo 

A nice feature of coercive functions is summarized in the fohowing 

THEOREM 4 (;i6l) Let f : W ^ TR. be continuous. If f is coercive, 
then f has at least one global minimizer. If, in addition, the first partial 
derivatives exist on IR"^, then f attains its global minimizers at its critical 
points. 

3 Unconstrained variational principles for the min- 
imal and maximal Br eigenvalues 

We now generahze the unconstrained variational principles of Auchmuty [2] 
to even order weakly symmetric tensors. We first consider the unconstrained 
optimization problem 

min/i(x) = —{Bx'^f + -Ax"^. (3.1) 

2m m 

When A and B are weakly symmetric, the gradient of the objective function 

/i is 

V/i(a;) = {Bx"')Bx'^'^ + Ax"'-^ . (3.2) 

The following theorem summarizes the properties of function /i . 

THEOREM 5 Assume that A and B are m -order n- dimensional real 
weakly symmetric tensors and B is positive definite. Let Amin be the smallest 
Br-eigenvalue of A. Then 

(a) /i is coercive on ]R°. 

(b) The critical points of f\ are 
(i) X = 0; and 

(a) any Br-eigenvector x of A associated with a Br-eigenvalue A < 0/ 



A satisfying Bx^ = — A. 

(c) If Amin < 0, then /i attains its global minimal value 

mm/i(x) = -7— A^in 
at any Br-eigenvector associated with the Br -eigenvalue Xmm satisfying Bx^ = 

(d) If Amin ^ 0, then x = is the unique critical point of /i and the unique 
global minimizer of fi on ]R°. 

Proof : (a) Since B is weakly symmetric positive definite, Theorem [3] 
asserts that 

^X™ > ;U||X||™, 

where ;U > is the smallest H-eigenvalue of B. This implies 

2m m 

as ll^ll — )• 00. Thus, /i is coercive on IR"^. 

(b) At a critical point of /i, its gradient V/i(x) = 0, that is, 

Ax"^-^ = -{Bx'^)Bx'^-\ (3.3) 

Clearly, x = is a critical point of /i as V/i(0) = 0. Moreover, if A < is 
a i3r.-eigenvalue of A, then 

Ax""-^ = XBx'^-K 

If X G ]R'^\{0} is a ^S^-eigenvector associated with this A and satisfies Bx"^ = 
—A, then it is a critical point of /i. 

(c) From Theorem [2l A > Amiru VA G cri3^{A). At the critical point x G 
IR"\{0} that is a ;S,.-eigenvector associated with a i3,.-eigenvalue A < and 
satisfies Bx"^ = —A, Ax"^ = — A^. Moreover, 

/i(^) = :^A2 - -A^ = -^A^ > -^ALn, 
Im m Zm 2m 

since > A > Amin- According to Theorem [4] and part (b), /i attains the 

1 2 

global minimum value Amin ^t any Sr-eigenvector associated with the 

2m 
5r. -eigenvalue Amin satisfying Bx"^ = — Amin- 

(d) Amin > implies that A > for any A G o"b^(^). Thus, Bx"^ = —A does 
not hold for any ;Br-eigenvector x of ^ associated with a ;Sr-eigenvalue A 



of A, as Bx'^ > for any x G ]R'^\{0} by the positive definiteness of B. 
Hence, x = is the unique critical point of /i. It is also the unique global 
minimizer of /i according to Theorem [H n 

We next consider the unconstrained optimization problem 

min/2(x) = —{Bx"'f - -Ax"^. (3.4) 

2m m 

When A and B are weakly symmetric, the gradient of the objective function 

/2 is 

Vf^{x) = {Bx"')Bx"'-^ - Ax"'-\ (3.5) 

Using a similar argument in the proof of the properties of /i, we can prove 
the following properties about f2- 

THEOREM 6 Assume that A and B are m -order n- dimensional real 
weakly symmetric tensors and B is positive definite. Let Amax ^e the largest 
Bj-- eigenvalue of A. Then 
(a) /2 is coercive on Ji^. 
(h) The critical points of /2 are at 
(i) X = 0; and 

(a) any Br-eigenvector x of A associated with a Br^eigenvalue X > of 
A satisfying Bx'^ = A. 

(c) If Amax > 0, then /2 attains its global minimal value 

min/2(x) = -^-^max 
at any B-eigenvector associated with the Br~eigenvalue Amax satisfying Bx"^ = 

^max- 

(d) //Amax < 0, then x = is the unique critical point 0//2. Moreover, it 
is the unique global minimizer 0//2 on TR^. 

Note that the functions /i and /2 are polynomials of degree 2m. A 
global polynomial optimization solver such as GloptiPoly3 [9j can be used 
to find the smallest Sf-eigenvalue of A and the largest iS^-eigenvalue of 
A by solving Problem (I3.ip and Problem (I3.4p respectively, provided that 
Amin < and Amax > 0. If Amin > 0, solving Problem ()3.ip does not result 
in the smallest Amin- In this case, we can solve the shifted problem 

min si(x, t) = — (i3x'")2 + —(A + tB)x'^, (3.6) 

xgR'i 2m m 



for a suitable t < 0. lit < — Aminj then the global minimum value of Problem 

p.6p is (Amin + t) . Thus, Amin Can be obtained by finding the global 

2m 
minimum of Problem ()3.6p . When Amax < 0, we can similarly solve the 

shifted problem 

min S2(x, t) = — (Sx™)2 -—(A + tB)x"', (3.7) 

xgR" 2m m 

for a suitable t > —Amax > to find Amax- 

Problem (jS.ip can be used to determine whether an even order symmetric 



tensor A is positive semidefinite or not. Take B = I ot B = In ■ If the 
global minimum value of /i equals 0, then A is positive semidefinite (or 
definite); otherwise, it is not. Assume that we have been able to determine 
that A is positive semidefinite. To further determine whether A is positive 
definite or semidefinite, we can solve ()3.6I) using t = —1. If the global 
minimum of si is —2m' then A is only positive semidefinite; otherwise A is 
positive definite. 

Local unconstrained optimization methods can be used to solve Problems 
p.ip and (j3.4p . These methods do not guarantee finding a global minimum. 
However, they converge to a critical point (see for example, [i5\). According 
to Theorems [5] and [6l the found nonzero critical point corresponds to a Br- 
eigenvalue of A. Therefore, local optimization solvers have the ability to find 
other eigenvalues besides the extreme ones. Moreover, if solving Problem 
(j3.ip with B = I^ or B = I results in a nonzero critical point, then it 



corresponds to a negative Z-eigenvalue or H-eigenvalue. This implies that 
local unconstrained optimization solvers can be used to solve Problem (|3.ip 
or Problem ()3.6p to determine if A is positive semidefinite. Finally, a local 
unconstrained optimization method such as the BFGS method has a fast 
rate of convergence - which is super linear. 

4 Numerical results 

In this section, we present some numerical results to illustrate the effec- 
tiveness of using the unconstrained variational principles for finding real 
eigenvalues of even order symmetric tensors. The experiments were done 
on a Laptop with an i3 CPU and a 4GB RAM, using MATLAB7.8.0 [H], 
the MATLAB Optimization Toolbox [14], and the Tensor Toolbox [3]. We 
did two groups of experiments: First, comparing the new approach with the 
SSHOPM method and the constrained optimization approach. Second, test- 



ing the ability of the new approach to determine positive semidefiniteness 
of even order symmetric tensors. 

4.1 Effectiveness of finding a Z-eigenvalue 

In our first group of experiments, we tested the new approach on finding 
Z-eigenvalues {B = I^ ) in order to compare it with the SSHOPM method 
(jllj). We will focus on solving Problem ()3.ip . The numerical behavior 
of solving Problem p.4p is similar. When B = I^ , the unconstrained 
variational principle (j3.ip becomes 



min/f (x) = —{x^xr + -Ax"". (4.1) 

2m m 

The gradient of the corresponding objective is 

V/f (x) = {x^x^-^x + Ax'^-K (4.2) 

We tested the symmetric 4 order tensors defined in the following examples: 

EXAMPLE 1 The 4 order n-dimensional tensor A is defined by 

-0.9, if i = j = k = l; 
A{i,j,k,l) = { l<i,j,k,l<n (4.3) 

0.1, otherwise. 

EXAMPLE 2 The 4*^^ order n-dimensional symmetric tensor A is gener- 
ated as follows: First randomly generate tensor T = randn(n, n,n, n), then 
use the symmetrize function in the Matlab Tensor Toolbox [3] to symmetrize 
T and obtain A = symmetrize(T). 

EXAMPLE 3 The 4 order n-dimensional symmetric tensor A is gen- 
erated as follows: First randomly generate tensor 3^ = randn(n, n,n,n); 
then create tensor Z by setting Z{i,j,k,l) = yTij^n', and finally use the 
symmetrize function in the Matlab Tensor Toolbox [3] to symmetrize Z and 
obtain A = symmetrize(2^). 

Since sometimes we are interested in finding the extreme eigenvalues of 
a tensor, we used the global polynomial optimization solver GloptiPolyS 
of Henrion, Lasserre, and Lofberg [9] to solve Problem (j4.ip for some 4 
order symmetric tensors in Examples 1-3. We observed that GloptiPolyS 
was able to solve (14. ip when n < 7. When n > 8, it was unable to solve 

10 



(|4.ip due to its memory requirement exceeding the capacity of the laptop 
computer we used. 

From now on in this subsection we shah focus on solving (|4.ip using 
a local optimization method. Specifically, we used the local optimization 
solver f minunc from the Matlab Optimization Toolbox |14j to solve Problem 
()4.ip . with its default settings except for the following: 

GradObj : on, LargeScale : off, TolX = TolFun = 10"^^, Maxlter = 1000. 

(4.4) 
We tested this approach and compared it with two other approaches on some 
tensors from Examples 1-3. 

4.1.1 Comparison with the constrained variational principle (jl.2p 

We tested and compared the unconstrained variational principle (j4.ip with 
the constrained variational approach (|1.2p (using B = I^ ) for finding Z- 
eigenvalues of some 4* order n-dimensional tensors given in Examples 1 and 
2. For the constrained variational principle approach, we used the fmincon 
function from the Matlab Optimization Toolbox [l3] to solve Problem ()1.2p , 
with its default settings except for the following: 

GradObj : on, GradConstr : on, TolX = TolFun = 10"^^, Maxlter = 1000. 

(4.5) 
Our numerical experiments indicated that Problem (jl.2p can be a sur- 
prisingly difficult problem for fmincon when ?7i = 4. Take the tensor A 
in Example [1] with n = 4 as an example. We ran f minunc on this tensor, 
using randomly generated initial vectors xq = randn(4, 1) 100 times and us- 
ing normalized randomly generated initial vectors xq = yo/\\yo\\2 100 times, 
where yo = randn(4, 1). We also ran fmincon on this tensor similarly. We 
observed that 



Solving (14. ip via f minunc: In all of 200 runs, this method successfully 
found the Z-eigenvalue A = —0.9345 of A. 

Solving (II. 2p with B = iZ via fmincon: (1) In the 100 runs using 
randomly generated initial vectors, this method successfully found the 
Z-eigenvalue A = —0.9345 of ^ 11 times and it failed to find a Z- 
eigenvalue of A in 89 runs. (2) In the 100 runs using normalized ran- 
domly generated initial vectors, it successfully found the Z-eigenvalue 
A = —0.9345 of A 43 times and it failed to find a Z-eigenvalue of A 
in 57 runs. The failures in both case (1) and case (2) were due to the 
divergence of the method. 



11 



4.1.2 Comparison with the SSHOPM method 

We now compare the unconstrained variational principle (|4.ip with the 
SSHOPM method of Kolda and Mayo [llj for finding Z-eigenvalues of some 
^th Qj-jjgj- n-dimensional tensors. The SSHOPM method is implemented as 
the sshopm function in the Matlab Tensor Toolbox p]. We used the default 
settings of sshopm except for the following: 

Tol = 10"^^, Maxlts = 5000. (4.6) 

The iterates x^ ' generated by sshopm keeps the norm ||x^ ' ||2 = 1- When 
f minunc converges to a nonzero critical point x of Problem ()4.ip correspond- 
ing to a Z-eigenvalue A < 0, ||5;||2^ = —A. To have a fair comparison of the 
two methods, we first normalized the nonzero vector x obtained by f minunc 
at termination so that x = x/||x||2. We define the error term by 

e=\\Ax^-Xx\\2, (4.7) 

where x is either the vector obtained by sshopm at termination or the nor- 
malized vector when solving (j4.ip via fminunc at termination. 

The SSHOPM method without shift (i.e., the original SHOPM method 
of Kofidis and Regalia |10] ) can fail to find a Z-eigenvalue, see for example, 
|10l 111] . We observed this behavior in our numerical experiments. Kolda 
and Mayo |llj proved that if the shift parameter a < is negative enough 
(or a > is large enough), then x'" generated by the SSHOPM method 
converges to a Z-eigenvector. However, the SSHOPM method slows down 
significantly when a very negative q < or very large a > is used [llj . 
Kolda and Mayo pT] found that using a = —2,— 1,1,2 worked well in their 
tests. 

Since the unconstrained variational principle (|4.ip leads to negative Z- 
eigenvalues of A, we used the SSHOPM method with a negative shift pa- 
rameter a. We solved ()4.ip via fminunc and ran sshopm with a = —2 on 
tensors of various dimensions from Example 1 and Example 3. For each 
tensor we tested, we ran each of fminunc and sshopm on the tensor 100 
times, using a normalized randomly generated initial vector 

^0 = 71 — fT' \^-^) 

\m\\2 

where yo = randn(n, 1) at each time. We report the numerical results in 
Tables 1 and 2, in which "CPU time" and "Accuracy" denote the "average 
CPU time (in seconds)" and "average e = \A^ — Ax||2" respectively. 



12 



Table 1: Unconstrained variational principle vs SSHOPM on some tensors 
from Example 1 using normalized randomly generated initial vectors 



Problem 


Method 


CPU time 


Accuracy 


n = 10 


SSHOPM {a = -2) 


0.16 


5.33 X 10- ' 


(14.111 via f minunc 


0.12 


2.12 X lO"'-" 


n = 20 


SSHOPM {a = -2) 


0.24 


4.64 X 10~' 


(|4.ip via f minunc 


0.18 


2.65 X 10"'' 


n = 30 


SSHOPM {a = -2) 


0.62 


3.73 X 10-' 


(14.111 via f minunc 


0.34 


3.70 X lO"'-" 


n = 40 


SSHOPM {a = -2) 


1.69 


4.05 X 10-' 


(14.111 via f minunc 


0.60 


3.10 X 10-" 


n = 50 


SSHOPM {a = -2) 


4.05 


3.63 X 10-' 


(14.111 via f minunc 


1.36 


3.30 X lO-'-" 


n = 60 


SSHOPM (a = -2) 


8.24 


3.08 X 10-' 


(|4.ip via f minunc 


2.69 


4.37 X lO-'-* 



Table 2: Unconstrained variational principle vs SSHOPM on some random 
tensors generated from Example 2 using normalized randomly generated 
initial vectors 



Problem 


Method 


CPU time 


Accuracy 


n = 10 


SSHOPM {a = -2) 


0.38 


1.03 X 10-" 


(14.111 via f minunc 


0.18 


9.98 X 10-" 


n = 20 


SSHOPM {a = -2) 


0.89 


1.33 X 10-" 


(|4.1|1 via f minunc 


0.27 


2.19 X 10-" 


n = 30 


SSHOPM {a = -2) 


2.31 


1.48 X 10-" 


(14.111 via f minunc 


0.53 


3.65 X 10-"* 


n = 40 


SSHOPM {a = -2) 


5.87 


1.61 X 10-" 


(|4.1|1 via f minunc 


1.06 


6.52 X 10-" 


n = 50 


SSHOPM {a = -2) 


11.84 


1.69 X 10-" 


(14.111 via f minunc 


2.21 


8.88 X 10-" 


n = 60 


SSHOPM {a = -2) 


24.28 


1.79 X 10-" 


(|4.ip via f minunc 


4.56 


1.09 X 10-' 
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Table 3: Unconstrained variational principle vs SSHOPM on a tensor gen- 
erated from Example 3 with n = 25, using normalized randomly generated 
initial vectors 



Method 


Min/Max/Mean Accuracy 


Min/Max/Mean CPU time 


(14.11) via fminunc 


1.06 X 10^73.44 X 10^74.23 x lO"-* 


0.36/1.14/0.63 


SSHOPM {a = 


.0) 


1.36 X 10^/4.02 X 10^/2.91 x 10" 


25.78/26.85/26.09 


SSHOPM (a = 


-1) 


1.08 X 10'*/5.26 X 10''/3.73 x 10" 


25.96/28.58/26.78 


SSHOPM (a = 


-2) 


2.56 X 10*/5.09 X 10*/3.97 x 10* 


25.84/27.50/26.34 


SSHOPM (a = 


-5) 


7.30 X 10V5.55 X 10^/3.58 x 10* 


25.79/27.11/26.25 


SSHOPM (a = 


-10) 


3.69 X 1074.55 X 10'*/1.56 X 10" 


25.78/25.91/25.81 


SSHOPM (a = 


-100) 


1.53/4.12 X 10*/2.80 x 10* 


25.66/ 25.84 /25.75 


SSHOPM (a = - 


-1000) 


6.50 X 10-71.09 X 10-78.79 X 10-' 


3.59/4.26/3.81 



Prom Tables 1 and 2, we observe that although both the SSHOPM 
method and the unconstrained optimization principle (j4.ip successfully find 
a Z-eigenvalue of A in all cases, solving (j4.ip via fminunc on average uses 
less CPU time than sshopm, particularly when n is large. This is perhaps 
due to the superlinear convergence property of the BFGS algorithm used in 
fminunc. 

A natural question is how to choose a suitable shift parameter a in 
sshopm. Using a = —2 worked well for the problems considered in Tables 
1 and 2. However, this choice is not a suitable one for some tensors from 
Example 3. We illustrate this in Table 3, in which we report the numerical 
results on a tensor from Example 3 with dimension n = 25. We ran sshopm 
with various shift parameters and fminunc on this tensor 10 times, using 
a normalized randomly generated initial vector xq as defined in (I4.8P with 
n = 25 at each time. The notations are: 

• Min/Max/Mean Accuracy denotes the minimum, maximum, and mean 

e = \\Ax^ — Xx\\2. 

• Min/Max/Mean CPU time denotes the minimum, maximum, and 
mean CPU time used. 

Clearly for this example, a = 0, —1, —2, —5, —10, —100 are not a suitable 
choice for the shift parameter. Although the tensor used in this test is arti- 
ficial, the numerical results indicate that choosing a suitable shift parameter 
can be crucial for the success of the SSHOPM method. 

We now summarize our comparison of the SSHOPM method and the 
unconstrained variational principle approach: 

• The SSHOPM method can be used to find a complex Z-eigenvalue and 
can handle both even and odd order symmetric tensors |11] . 
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• The unconstrained variational principle approach (e.g. solving (j4.ip 
via fminunc) can be faster than the SSHOPM method particularly 
when n is large. 

• If the purpose is to find one Z-eigenvalue for a given tensor, choosing 
a suitable shift parameter can be crucial for the SSHOPM method. 
Although sometimes it needs to solve a shifted problem (|3.6|) or (j3.7|) 
(see the next subsection), choosing a suitable shift parameter is less 
critical for the unconstrained variational principle approach. 

• A global polynomial optimization solver can be used to solve the op- 
timization problems arisen in the unconstrained variational principles 
to find the largest or smallest eigenvalues. The found Z-eigenvalue by 
the SSHOPM method is not necessarily the largest or the smallest one. 

4.2 Determining positive semidefiniteness 

In some applications, it is important to determine if an even order symmetric 
tensor A is positive semidefinite (see, for example, [211 EZ])- An attractive 
property of function ()3.ip is that any of its nonzero critical point is a Br- 
eigenvector of A corresponding to a ;Sr.-eigenvalue A < 0. This feature allows 
us to use a local optimization solver to determine the positive semidefinite- 
ness of an even order symmetric tensor, since A is positive semidefinite if 
all of its H-eigenvalues and all of its Z-eigenvalues are nonnegative ( |17j ) . 

To test the effectiveness of using the unconstrained optimization ap- 
proach to determine the positive semidefiniteness of even order symmetric 
tensors, we did some preliminary numerical experiments in which Problem 
()3.6p was solved via fminunc. The parameters for fminunc were the same 
as in g3D. Note that ([M]) becomes (^ when t = 0. When t / 0, solv- 
ing ()3.6p leads to a ^Br-eigenvalue of the tensor A + tB. In this situation, 
subtracting t from the computed eigenvalue oi A + tB will result in a Br- 
eigenvalue of A. We tested both B = ^(corresponding to H-eigenvalues) 
and B = /™ (corresponding to Z-eigenvalues). We summarize the uncon- 
strained optimization approach for determining the positive semidefiniteness 
of an even order symmetric tensor A in Algorithm 1. 
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ALGORITHM 1 

Input: Tensor A. 

Step 0. Choose parameters t < 0, 772 > ^i > 0, and tensor B = X or 

B = C'\ 

Step 1. Solve Problem (|3.6p with parameter t approximately using 

a local optimization solver such as fminunc. Let x and s denote the 

approximate optimal solution and optimal objective value of (13. 6p found 

by the optimization solver. 

Step 2. 

• If Bx'^ < rji, then A is positive definite, stop. 



If Bx"^ > r/2, then set A = —\/—2ms — t (which is a i3j.-eigenvalue 
of ^). If A > 0, then A is positive semidefinite; otherwise A is not 
positive semidefinite, stop. 



REMARK 2 If ??i < Bx^ < r]2, then Algorithm 1 is inconclusive. In this 
case, we can use a different shift parameter t < and repeat Step 1. 

For comparison, we also tested the SSHOPM method. For this method, 
we used the same setting as in ()4.6p except for that Maxlts was changed to 
Maxlts = 10000. 

EXAMPLE 4 The 4*^ order n-dimensional symmetric tensor A is gener- 
ated as follows: First randomly generate tensor T = randn(n, n,n,n), then 
use the symmetrize function in the Matlab Tensor Toolbox [3J to symmetrize 
T and obtain Z = symmetrize(7'). Finally set 

1000, ifl<i = j = A; = /<n-l; 

A{i,j,k,l) = { -1, [ii = j = k = l = n; (4.9) 

Z{i,j,k,l), otherwise. 

A is not positive semidefinite when n > 2. 

EXAMPLE 5 The 4*^^ order 3-dimensional tensor A is defined by ^(1, 1, 1, 1) 
1; ^(2,2,2,2) = 0; ^(3,3,3,3) = -0.001; and A{i,j,k,l) = for ah other 
i,j, k, I. This tensor is not positive semidefinite. 
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We tested Algorithm 1 using B = I or B = I™ and shift parameters 
t = or t = — 1 and compared them with the SSHOPM method with 
different shift parameters on some tensors from Examples 4 and 5. For 
each tensor, we ran each of these methods 100 times, using a normalized 
randomly generated initial vector as defined in (|4.8p at each time. We used 
rji = 10^^^ and rj2 = 10~^ in Algorithm 1. 

We report the numerical results in Tables 4 and 5. In both tables, 
"Success rate" denotes the percentage of times where a negative eigenvalue 
was found (and therefore the corresponding method successfully determined 
that A is not positive semidefinite) ; "CPU time" denotes the average CPU 
time (in seconds); and "NIT" denotes the average number of iterations used 
by the SSHOPM method. 

For the tensor generated from Example 4, the SSHOPM method always 
converged to a positive eigenvalue when a = —2, a = — 10, a = —50. When 
a = —100, it converged to a negative eigenvalue in 34 out of 100 runs. 
It successfully found a negative eigenvalue when a = —500 in all of the 
100 runs, using an average CPU time of 30 seconds. On the other hand. 
Algorithm 1 using B = Z or B = I^ and t = Oort = — 1 successfully 
found a negative eigenvalue in all runs, using much less CPU time. 

For the tensor from Example 5, when t = 0, Algorithm 1 using B = X cor- 
rectly identified that A is not positive semidefinite 91% of times; Algorithm 1 
using B = In was only successful 32% of times. The failures in both cases 
were due to that fminunc converged to the critical number re = of (j3.ip . 
This is because fminunc is a local optimization solver. It only guarantees 
to converge to a critical point. We found that using a negative parameter t 
can significantly increase the success rate, particularly in the B = /^ case. 
Indeed, when t = —1 was used, Algorithm 1 successfully found a negative 
eigenvalue 98% of times in both B = X and B = I^ cases. The two failure 
runs in each case were due to that the eigenvalue A = was found. In all 100 
runs, the SSHOPM method (with a = —2) terminated when the maximum 
allowed number of iterations (which is 10000) was reached. In 10 out of 100 
runs, the method terminated at an approximate Z-eigenvalue very close to 
A = (and hence failed to correctly determine the positive semidefiniteness 
of ^). We plot the computed Z-eigenvalues by the SSHOPM method in the 
100 runs in Figure 1. From this figure we observe that the SSHOPM method 
successfully determined that A is not positive semidefinite in less than 90% 
of times. 

In summary. Algorithm 1 using B = I or B = In and a negative shift 
parameter t is more efficient than the SSHOPM method on determining the 
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Table 4: Determining positive semidefiniteness using a tensor generated in 
Example HI n = 30 



Method 


Shift parameter 


Success rate 


CPU time 


NIT 


Algorithm 1 (B = I) 


t = 


100 


1.36 




Algorithm 1 (B = X) 


t = -1 


100 


1.26 


Algorithm 1 (B = C'") 


f = 


100 


0.99 




Algorithm 1 (B = /,?'") 


t = -1 


100 


0.99 


SSHOPM 


a = -2 





2.97 


409.4 


SSHOPM 


a = -10 





2.34 


310.1 


SSHOPM 


a = -50 





7.41 


1111.2 


SSHOPM 


a = -100 


34 


3.22 


449.3 


SSHOPM 


a = -500 


100 


30.00 


4693.4 



Table 5: Determining positive semidefiniteness using the tensor in Example 



Method 


Shift parameter 


Success rate 


CPU time 


NIT 


Algorithm 1 (B = X) 


i = 


91 


0.28 




Algorithm 1 (6 = 1) 


t= -1 


98 


0.62 


Algorithm 1 (B = C''^) 


t = 


32 


0.16 




Algorithm 1 (B = ir'") 


t = -1 


98 


0.47 


SSHOPM 


a = -2 


<90 


34.58 


10000 



positive semidefiniteness of the tensors from Examples 4 and 5 we tested in 
terms of CPU time. 

We now turn to the performance of Algorithm 1 on positive semidefinite 
(definite) tensors. We tested some tensors from Examples [6] and [7] using 
Algorithm 1 with B = T and t = — 1 and with B = I^ and t = — 1. For 
each tensor we tested, we run each method 100 times, using a normalized 
randomly generated initial vector as defined in ()4.8p at each time. We report 
the numerical results in Tables 6 and 7. In Table 6, "min Bx"^" denotes the 
smallest Bx"^ in 100 runs and "Success rate" denotes the percentage of times 
where the minimum eigenvalue A = was obtained. In Table 7, "min/max 
i3x™" denotes the smallest and largest Bx"^ in 100 runs and "Success rate" 
denotes the percentage of times when the method correctly determined that 
A is positive definite. From Tables 6 and 7, we observe that Algorithm 1 
using ";B = I and t = —1" or ";S = I^ and t = —1" was able to efficiently 
determine the positive semidefiniteness of the tensors from Examples 6 and 
7 we tested. 



EXAMPLE 6 The 4 order n-dimensional tensor A is defined by A{k, k, k, k) 
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Figure 1: The computed Z-eigenvalues by the SSHOPM method in the 100 
runs on tensor A from Example 5. 



rand(l) for k = 1,2, • • • ,n — 1; A{n,n,n,n) = 0; and A{i,j, k,l) =0 for all 
other i,j, k, I. A is positive semidefinite. 

EXAMPLE 7 Consider the positive definite 4*^^ order n-dimensional ten- 
sor A defined by A{k, k, k, k) = 10k for A; = 1, 2, • • • , n; and A{i,j, k,l) = 
for all other i,j, k,l. 



Table 6: Determining positive semidefiniteness using a tensor in Example [U 
n = 30 



Method 


mill Bx"" 


Success rate 


CPU time 


Algorithm 1 {B = I,t = - 


-1) 


1.00 


100 


1.68 


Algorithm 1 (B = C^', t = 


-1) 


1.00 


100 


2.54 



Table 7: Determining positive semidefiniteness using a tensor in Example [71 
ra = 30 



Method 


min/max Si™ 


Success rate 


CPU time 


Algorithm 1 {B = I,t = -1) 


7.12 X lO"'' / 2.65 X 10"'" 


100 


1.66 


Algorithm 1 (B = C^', t = -1) 


4.98 X 10"'^ / 1.97 X 10"^^ 


100 


1.81 



19 



5 Final Remarks 

We have introduced two unconstrained optimization problems and obtained 
some variational characterizations for the minimum and maximum Br eigen- 
values of an even order weakly symmetric tensor, where B is weakly sym- 
metric positive definite. These unconstrained optimization problems can 
be solved using some powerful optimization algorithms, such as the BFGS 
method. We have provided some numerical results indicating that our ap- 
proach of solving Problem (|4.ip via f minunc compares favorably to the ap- 
proach of solving (|1.2|) via f mincon and the SSHOPM method for finding a 
Z-eigenvalue of an even order symmetric tensor. This approach can also be 
used to find other Br eigenvalues of even order symmetric tensors, includ- 
ing H-eigenvalues and D-eigenvalues and to find a Br eigenvalue of an even 
order weakly symmetric tensor. Furthermore, we have provided some nu- 
merical results that show this approach is promising on determining positive 
semidefiniteness of an even order symmetric tensor. 

A direction for future research is to develop a global polynomial opti- 
mization algorithm that can solve problems (|3.ip and (j3.4p with larger m 
and n. 
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